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Abstract 

This paper presents the large-eddy simulation of the lid-driven cubic cavity flow by 
the spectral element method (SEM) using the dynamic model. Two spectral filtering 
techniques suitable for these simulations have been implemented. Numerical results 
for Reynolds number Re = 12'000 are showing very good agreement with other 
experimental and DNS results found in the literature. 
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1 Introduction 



Spectral element methods have been mainly applied to the direct numerical 
simulation (DNS) of fluid flow problems at low and moderate Reynolds (Re) 
numbers. With the advent of more powerful computers, especially through 
cluster technology, higher Re values seem to fall within the realm of feasibility. 
However, despite their high accuracy, spectral element methods are still far 
from reaching industrial applications that involve developed turbulence at Re 
values of the order of 10^ — 10^. The reason for that dismal performance is that 
a resolved DNS including all scales from the gross structures to Kolmogorov 
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scales, needs a number of degrees of freedom (dof) that grows like Re^^^. 
Therefore with increasing Re, we have to increase the number of elements, E, 
and the degree, N, of the polynomial spaces. This places the computational 
load far out of the reach of present day computers. 

Large-eddy simulation (LES) represents an alternative to DNS insofar that it 
involves less dofs because the behaviour of the small scales are modeled. The 
LES methodology deals with high Re number unsteady flows by using coarse 
meshes in which the contributions of small sub-grid scale (SGS) structures are 
modeled while the large-scale structures are obtained by the computed flow 
dynamics. The reader is referred to the monograph by Sagaut [15] for details. 
In this paper we will focus our attention on the dynamic model [4, 11, 15] 
which will be evaluated by comparing LES numerical results and DNS on the 
three-dimensional cubic cavity flow. 



2 Mathematical modelling 

2.1 The governing equations 

Large-scale quantities, designated by an "overbar" , are obtained by a filtering 
process on the domain Q. Assuming the filter commutes with differentiation 
and applying the filter to the Navier-Stokes equations in divergence form for 
the nonlinear term, one obtains the following relations: 



— V- vv = -Vp-Fi/Av- V-T, (1) 

at 

V-v = 0. (2) 

Here, v is the filtered velocity, t denotes the time, p is the filtered pressure 
divided by the constant density, u the kinematic viscosity. The symbols V and 
A represent the nabla and Laplacian operators, respectively. The SGS stress 
tensor r takes the small-scale effects into account and is given by 

T = VV — V V. (3) 

2.2 Smagorinksy and dynamic models 

The SGS Smagorinsky model (SM) [17] uses the concept of turbulent viscos- 
ity and assumes that the small scales are in equilibrium, balancing energy 
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production and dissipation. This yields the following expression for the eddy- 
viscosity 

UT = {CsAf\S\, (4) 

where |S| = {2SijSijY/'^ is the magnitude filtered strain-rate-tensor, Cs is the 
Smagorinsky constant and A the filter width. The SM has several drawbacks. 
The most severe one is the constant value of Cs during the computation which 
produces too much dissipation. Furthermore the SM does not provide the 
modeller with backscattering where kinetic energy is transferred from small 
scales to larger scales. 

The dynamic model (DM) proposed by Germano et al. [4] overcomes the 
difficulty of constant Cs, by allowing it to become dependent of space and 
time. Now we have a dynamic constant = Cd(x, t). Let us introduce a test- 
filter length scale A that is larger than the grid length scale A (e.g. A = 2A). 
Using the information provided by those two filters and assuming that in the 
inertial range of the turbulence energy spectrum, the statistical self-similarity 
applies, we can better determine the features of the SGS stress. With the test 
filter, the former LES equations ([1]) yield a relation involving the sub-test-scale 
stress 

T = W-VV. (5) 

We introduce the Germano identity to obtain the relation between T and the 
filtered r such that 

L = T-f = VV-VV. (6) 

We apply the eddy viscosity model to r and T and we obtain using the self- 
similarity hypothesis for the constant Ca 

T - ^tr(T) I = -2CdA'|S|S = Cd/3, (7) 

and 

T - ^tr (T) I = -2Cd^^ 1 1| t = Cdtt, (8) 

where the symbol tr denotes the trace of the tensor. Inserting ([7]) and (IE]) in 
the deviatoric part L'^ of L produces 

L - ^tr(L)I = = CdCK - 6^(3. (9) 

Assuming that Cd does not vary too much in space, one sets ~ Cd and 
one can deduce from a least square minimization of the error related to ([9]) 
(see [12,15]) that 

= '°r^>-'-^. (10) 

(a-/3):(o-/3)' 

where the notation : is used for inner tensor product (double contraction). 
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3 Numerical approximation 



3.1 Space discretization 

The numerical approximation is obtained through a weak formulation of the 
Eqs. ([I])-© discretized using the Lagrange-Legendre approximation. The 
reader is referred to the monograph by Deville, Fischer and Mund [3] for de- 
tails. The velocity and pressure are expressed in the Pat — PAr-2 spaces where 
Pat is the set of polynomials of degree < in each space direction. This ap- 
proximation avoids the presence of spurious pressure modes as it was proved 
by Maday and Patera [13]. The quadrature rules involved in the weak for- 
mulation define a Gauss-Lobatto-Legendre (GLL) grid for the velocity nodes 
and a Gauss-Legendre grid (GL) for the nodal pressures. 

Borrowing the notation from [3] , the semi-discrete filtered Navier-Stokes equa- 
tions resulting from the spectral element discretization are 



The diagonal mass matrix M is composed of d blocks, the mass matrices M, 
with (i = 3 for the three-dimensional cavity problem. The supervector v con- 
tains all the nodal velocity components while p is made of all nodal pressures. 
The matrices K, D^, D are the discrete Laplacian, gradient and divergence op- 
erators, respectively. The matrix operator C represents the action of the non 
linear advection term written in convective form v ■ V, on the velocity field 
and depends on v itself. The spatial discretization leads to a set of non linear 
ordinary differential equations ffTTl) subject to the incompressibility condition 



3.2 Time integration 

As the LES viscosity is not invariant, we modify the standard time integration 
scheme in such a way that this space varying viscosity be handled explicitly 
as this was done e.g. in [1,6]. Let us define the effective viscosity as 



where Ucst is the sum of v and the spatial average of ut over the computational 
domain, being by construction constant in space but not in time. The filtered 
Navier-Stokes equations become 



M + C V + z/ K V - D^p + Dr = 0, 




(11) 
(12) 



(13) 
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M-^ + z/est Kv - D^p = -Cv + 2D(i.efr - z/cst)S (14) 
-Dv = 0. (15) 

The viscous linear term and the pressure are imphcitly integrated by a back- 
ward differentiation formula of order 2 (BDF2) while the terms in the right- 
hand side of Eq. f|T^ are computed by a second order extrapolation method 
(EX2). The explicit viscous part 2D(z/efr — z/cst)S leads to a stability condition 
such that (z/cff — z/cst)^^ ^ C /N^ while the CFL condition restricts the time 
step such that u^a^/S.t < C/N"^. It would seem that the viscous restriction 
is more severe than the convective one. However the magnitude of the term 
2D(z/eff — 1^081)3. is far smaller than the one of the convective term. Therefore 
the stability is indeed enforced by the CFL limit. 

The implicit part is solved by a generalized block LU decomposition, using a 
standard fractional-step method with pressure correction which may be pre- 
conditioned by various algorithms. 



4 Filtering 

As spectral elements offer high accuracy for the flow at hand, we construct 
the filters using two spectral techniques. The first one is a nodal filter acting 
in physical space on the nodal velocity components (and pressure) to stabilize 
the computations. The second method is designed as a modal filter and is 
carried out element-wise in spectral space and corresponds to the convolution 
kernel of the LES filtering. 

4-1 Nodal filter 

The nodal filter is due to Mullen and Fischer [14] and is adequately suited 
to parallel spectral element computation. Introducing hNj,j = 0,. . . ,N the 
set of Lagrange-Legendre interpolant polynomials of degree on the GLL 
grid nodes ^N,k, k = 0, . . . ,N, the rectangular matrix operator /^^ of size 
(M + 1) X (A^ + 1) is such that 

Therefore, the matrix operator of order A^ — 1 

n^v-i = i^-' (17) 

interpolates on the GLL grid of degree A^— 1 a function defined on the GLL grid 
of degree A^ and transfers it back to the original grid. By this process, one can 
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show that one gets rid of the highest modes of the polynomial representation. 
The ID filter is given by the relation 

v=[aUN-i + {l~a)I^]v. (18) 

The LES version of the filter sets a = 1 and is given by 

v = ImInV, (19) 

where M is equal to — 2 or iV — 3. The extension to three-dimensional 
problems results easily from the matrix tensor product properties of the filter. 



4.2 Modal filter 



The variable v aside its Lagrange- Legendre representation may also be ap- 
proximated by a modal basis that was first proposed in the p version of the 
finite elements. That basis was used by Boyd [2] as a filter technique and is 
built up on the reference parent element as 

00 = 4>i = 0fc = LkiO - Lk-2{0, 2<k<N. (20) 

The one-to-one correspondence between the Lagrangian basis and the p rep- 
resentation yields 

N 

vi^^)=T.^kM^^) (21) 
k=0 

which in matrix notation reads 

V = $v. (22) 

The filter operation is performed in spectral space through a diagonal matrix 
T with components such that 

= (1 + il/Ky) ^^^^ 

where the cut-off value kc corresponds to = 1/2. The entire filtering process 
for a. ID problem is given by 

v = G'^v = $T$"^v (24) 

The extension to three-dimensional is trivial by the matrix tensor product 
properties. 

The transfer function of Eq. (123|) corresponds to a standard low-pass filter 
required when using any spectral method [2]. 
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4-3 The filter length 



For a ID problem using the spectral element method, a common choice [7] of 
filter length is 

A = (25) 
P 

where s is the element size and p the highest polynomial degree in the spectral 
decomposition Eq. fl2Tl) that is the closest to the frequency kc 

p = k, such that inf(| k-k^l) <l/2, k = 0,...,N. (26) 

We notice that the filter length decreases when the element is refined. The 3D 
formula for rectilinear elements is 

A{x,y,z) = {A^{x)A,{y)A,{z)Y/' = (^^'-l] ' . (27) 

\P\ Pi P3 J 



5 The lid-driven cubical cavity problem 



The lid-driven cavity presents although the geometry is simple, complex phys- 
ical phenomena. As in this case, we have no homogeneous direction, the pres- 
ence of side walls confining the full flow modifies the flow patterns and the 
route to turbulence. For the description of the physics among the abundant 
literature, we refer the reader to Koseff and Street [8], Jordan and Ragab [5], 
Leriche and Gavrilakis [10], and Shankar and Deshpande [16]. 

Figure [Hb shows the cubical cavity. The flow motion is induced by the top 
lid that moves in the x-direction with a constant unit velocity Uq = 1. The 
Reynolds number is consequently Re = f/o2/i/z/. We will essentially address 
the case of the flow at Re = 12'000. The kinetic energy is provided to the flow 
by the shear stress at the top lid through viscous diffusion. The amplitude of 
the Reynolds stress below the lid is negligible indicating that the flow under 
the lid is mainly laminar but transient. The momentum transfer from the lid 
induces a region of strong pressure in the upper corner of the downstream 
wall as the flow, mainly horizontal prior the corner, has to change direction 
and moves vertically downwards. This sharp turn dissipates energy in that 
region. Along the downstream wall the plunging flow behaves like a jet with 
a variable thickness. Near the symmetry plane the jet thickness is reduced 
while it increases away from this plane. This jet, laminar and unsteady at 
the very beginning, separates from the cavity wall at mid-height and grows as 
two elliptical jets on both sides of the symmetry plane. They hit the bottom 
wall where they produce turbulence. This turbulence is convected away by the 
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main central vortex towards the upstream wall where the flow slows down and 
relaminarizes during the fluid rise. 



>■ 0- 



-0.5- 



■'-1 its i5 M 

" -h h " 

Fig. 1. (a) (left) Spectral element grid in the mid-plane z = and (6) (right) Sketch 
of the geometry of the lid-driven cubical cavity. 

In order to resolve the boundary layers along the lid and the downstream wall, 
the spectral elements are unevenly distributed as can be seen in Figure [Ha. 
The spatial discretization has = Ey = = 8 elements in the three space 
directions with = Ny = = 8 polynomial degree. The spectral element 
calculation has two times less points per space direction than the DNS of 
Leriche-Gavrilakis [10] who employed a 129'^ Chebyshev discretization. 

As far as the velocity imposed on the lid is concerned, the unit velocity induces 
severe discontinuities along the top edges. In order to remove these defects we 
use a high degree polynomial as Leriche and Gavrilakis did 

,v = w = 0. (28) 
No-slip conditions are applied to the other walls. 

Both nodal and modal fllters were used in our LES computations; the former 
with M = — 2 to stabilize the velocity fleld at each time step and the 
latter with = N — 2 to fllter the highest modes in the Legendre space. The 
computations are particularly sensitive to the values of M and k^, smaller 
values will affect spectral convergence whereas higher values will have very 
little effect on the smallest scales of the problem. 

The dynamic constant produced has high values in the regions of high 
velocity gradients. Its maximum value fluctuates around 0.25 with locally some 
negative values which are eliminated by clipping. 

The time step is chosen as At = 10~^ and the complete simulation comprises 
201'000 iterations leading to a total effective simulation time of 201 time units. 
The reference results are the DNS data of Leriche [9] and the experimental 
ones from Koseff and Street [8], corresponding to I'OOO and 145.5 time units 
respectively. In the cavity flow, the average is obtained by time averaging. 




u 



1 - {x/hf^ ' 1 - {z/h) 



,18 



158 



The results of an under-resolved DNS performed on the coarse grid made of 65^ 
grid points and designed for the LES are presented on Fig. [2lb. These results 
are compared to the same results for the LES (Fig. [2Ja) and for the DNS (Fig. 
IHc). A qualitative comparisons of the three sets of mean- velocity contours for 
(f/) and {V) allows to conclude to the failure of the under-resolved DNS to 
reproduce the physics involved in this problem. As expected the resolved-scale 
spectrum is smaller than the full-scale spectrum for this turbulent flow and 
therefore necessitates to model the under-resolved scales as it is done in our 
LES (Fig. Ha). 




Fig. 2. Contours of the x-component of the mean velocity field (U) (row (i)) and 
of the y-component {V) (row (ii)) in the mid-plane z = for column (a) the 
LES, column (6) the under-resolved DNS, column (c) the DNS from Leriche and 
Gavrilakis [10]. 



From a more quantitative point-of-view, the x-component u' and the y-compo- 
nent v' of the rms fluctuations of the velocity field have been computed along 
two ID lines x = and y = in the mid-plane z = 0. Results for the LES are 
shown on Figure [3] together with the DNS and experimental results, proving 
the excellent performance of our SGS modelling. The differences of the three 
sets of data in terms of peak amplitude observed can be explained by the three 
different ensemble averaging used as proved by Leriche [9] . A larger ensemble 
averaging tends to reduce the amplitude of the different peaks. 
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(a) rms fluctuations u' jJJ 



(6) rms fluctuations u' /U 




Fig. 3. Profiles of the rms fluctuations u' ((a) and (6)), v' ((c) and (d)) in the 
mid-plane z = 0, along the lines x = ((6) and {d)) and y = ((a) and (c)), for 
the LES, the DNS from Leriche and Gavrilakis [10] and the experimental data from 
Koseff and Street [8]. 

6 Conclusion and future studies 



A large-eddy simulation of the three-dimensional lid-driven cubical cavity flow 
using the spectral element method has been presented. The treatment of the 
subgrid-scales reUes on a dynamic model for the eddy viscosity. This LES has 
been carried out on a parallel architecture with a relatively coarse grid and 
numerical results appeared to be extremely close to the DNS and experimental 
results available in the literature. Moreover the results of the under-resolved 
DNS on the coarse grid are far from providing any insight into the physics of 
the problem. 

Our next goal is to perform the same LES but with other treatments of the 
subgrid-scales and determine the most efficient model when the simulation are 
based on the spectral element method. 
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